data{
  int<lower=1> P;//stimulation
  int<lower=1> J;//item
  int<lower=1> I;//respondent
  int<lower=1> K;//dim
  real delta[P,J,I];//data
  real obs_data[P,J,I];
}

parameters{
  matrix[J-3,K] Z_raw;
  matrix[P,K] X;
  vector<lower=0>[K] sigx;
	real<lower=0> sigd;
	vector<lower=0>[K] sigz;
	real<lower=0> p21;
	real p31;
	real<lower=0> p32;
}

transformed parameters{
  matrix[J,K] Z;
	for(j in 1:J){
	  for(k in 1:K){
	    if(j == 1 && k == 1){
	      Z[j,k] = 0;
	    }else if(j == 1 && k == 2){
	      Z[j,k] = 0;
	    }else if(j == 2 && k == 1){
	      Z[j,k] = p21;
	    }else if(j == 2 && k == 2){
	      Z[j,k] = 0;
	    }else if(j == 3 && k == 1){
	      Z[j,k] = p31;
	    }else if(j == 3 && k == 2){
	      Z[j,k] = p32;
	    }else{
	      Z[j,k] = Z_raw[j-3,k];
	    }
	  }
	}
}

model{
  //model
  for(i in 1:I){
    for(j in 1:J){
      for(p in 1:P){
        real temp = 0;
        if(obs_data[p,j,i] == 1){
          for(k in 1:K){
            temp += (X[p,k]-Z[j,k])^2;
          }
          delta[p,j,i] ~ normal(sqrt(temp),sigd);
        }
      }
    }
  }
  //coordinates
  for(k in 1:K) Z_raw[,k] ~ normal(0,sigz[k]);
  for(k in 1:K) X[,k] ~ normal(0,sigx[k]);
  //constrained coordinates
	p21 ~ normal(0,sigz[1]);
	p31 ~ normal(0,sigz[1]);
	p32 ~ normal(0, sigz[2]);
	//prior
  sigx ~ cauchy(0,5);
	sigd ~ cauchy(0,5);
	sigz ~ cauchy(0,5);
}


